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ABSTRACT 

Under the assumption of thermal wind balance and effective entropy mixing in 
constant rotation surfaces, the isorotational contours of the solar convective zone 
may be reproduced with great fidelity. Even at this early stage of development, 
this helioseismology fit may be used to put a lower bound on the midlatitude 
radial solar entropy gradient, which in good accord with standard mixing length 
theory. In this paper, we generalize this solar calculation to fully convective 
stars (and potentially planets), retaining the assumptions of thermal wind bal- 
ance and effective entropy mixing in isorotational surfaces. It is found that each 
isorotation contour is of the form R 2 = A + B$>(r), where R is the radius from 
the rotation axis, $(r) is the (assumed spherical) gravitational potential, and A 
and B are constant along the contour. This result is applied to simple models 
of fully convective stars. Both solar-like surface rotation profiles (angular veloc- 
ity decreasing toward the poles) as well as "antisolar" profiles (angular velocity 
increasing toward the poles) are modeled; the latter bear some suggestive resem- 
blance to numerical simulations. We also perform exploratory studies of zonal 
surface flows similar to those seen in Jupiter and Saturn. In addition to providing 
a practical framework for understanding the results of large scale numerical sim- 
ulations, our findings may also prove useful in dynamical calculations for which a 
simple but viable model for the background rotation profile in a convecting fluid 
is needed. Finally, our work bears directly on an important goal of the CoRoT 
program: to elucidate the internal structure of rotating, convecting stars. 

Subject headings: convection — hydrodynamics — stars: rotation — Sun: rota- 
tion — Sun: helioseismology 



1 Laboratoire de Radioastronomie, Ecole Normale Superieure, 24 rue Lhomond, 75231 Paris CEDEX 05, 
France steven.balbus@lra.ens.fr 

2 Adjunct Professor, Department of Astronomy, University of Virginia, Charlottesville VA 22903, USA 

3 DAMTP, Centre for Mathematical Sciences, Wilberforce Road, Cambridge CB3 OWA, UK 



-2 - 



1. Introduction 

Recent work (Balbus 2009, hereafter B09; Balbus et al. 2009, hereafter BBLW) suggests 
that, away from inner and outer boundary layers, the isorotation contours in the bulk of the 
solar convective zone (hereafter SCZ) correspond to the characteristic curves of the vorticity 
equation in its quasi-linear "thermal wind" form. Since both the entropy 5* and angular 
velocity Q figure in this single equation, the mathematical existence of these characteristics 
requires that there be some functional relationship between S and Q. 

B09 put forth the possibility that this relationship was a consequence of the SCZ being 
marginally stable to axisymmetric magnetobaroclinic modes. In BBLW, however, it was 
shown that an entropy/angular velocity relationship is expected to be present under very 
general conditions; a magnetic field is not a prerequisite. Because robust convective struc- 
tures lie in surfaces in which Q is constant, and because these structures mix entropy very 
efficiently, surfaces of constant entropy and surfaces of angular velocity correspond with one 
another. More precisely, it is not quite the entropy that is homogenized by the convective 
mixing, it is the "residual entropy:" some radial entropy profile must be present to trigger 
the convection, and it is the residual entropy, the entropy that departs from the underlying 
radial profile, that is homogenized by the convective mixing in constant Q surfaces (BBLW). 
The confluence of constant angular velocity and constant residual entropy surfaces ensures 
that there is a functional relationship between them, and the existence of this functional 
relationship, even in the absence of knowledge of its precise form, allows one to deduce the 
structure of the isorotation characteristics of the thermal wind equation (hereafter TWE). 
If isorotational surfaces are the same as constant residual entropy surfaces, this is a result of 
considerable practical interest: in Appendix II, we show that it may be used in conjunction 
with helioseismology data to place a lower bound on the logarithmic radial entropy gradient 
(an otherwise observationally inaccessible quantity) of a few times 10~ 6 . This number is 
good accord with mixing length theory estimates (e.g. Stix 2004)0. 

Because the calculations presented in this paper depend strongly on the assumptions of 
thermal wind balance and the confluence of constant rotation and residual entropy surfaces, 
we review the supportive arguments that apply in environments similar to the SCZ. There 
are four distinct strands: 

Numerical Simulations. The first piece of evidence is provided by numerical simulations. 
A comparison of figure (2e) (isorotational contours) and (3a) (residual entropy) in Miesch, 



*A recent helioseismic inversion calculation by Brun et al. (2009) reports very large entropy variations 
in the SCZ, inconsistent with previous numerical simulations, thermal wind balance, and standard mixing 
length theory. These results may, however, be noise dominated. 
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Brun, & Toomre (2006; see also figure [2] of BBLW) shows a similarity that is striking to the 
eye. By way of contrast, there is no such correspondence between surfaces of total entropy 
and angular velocity (e.g., figure [2] in BBLW). The validity of TWB is evident in figure (5) 
of Miesch et al. (2006). 

Dynamics. The second line of argument is based on a very simple physical picture put forth 
in BBLW, and has already been touched upon above. Localized, embedded disturbances in a 
shearing system will be distorted into surfaces of constant Q. This follows from mass conser- 
vation alone, and is most easily shown by using coordinates comoving with the background 
flow (BBLW). But the embedded disturbances that are of interest here are convective rolls, 
which mix (and thus homogenize) residual entropy. Since the convective rolls are at once in 
surfaces of constant rotation and constant residual entropy, these two types of surfaces must 
correspond. 

Agreement between theory and observations. Even if a functional relationship between the 
angular velocity and residual entropy were simply postulated ad hoc, the resulting agreement 
between the calculated isorotational contours and the helioseismology data is so striking that 
the fit itself is an independent piece of evidence for the reasoning on which it is based (cf. 
fig. [1] of BBLW). The relationship traansforms the TWE from an equation relating certain 
partial derivatives of Q and S into a self-contained, quasi-linear partial differential equation 
for Q alone. The isorotation contours then emerge as the characteristics of this equation. 
These characteristics thus depend for their very existence on a functional relationship be- 
tween (residual) entropy and angular velocity. Without a functional relationship between 
entropy and angular velocity, the notion of isorotation characteristics would make no math- 
ematical sense. With the functional relationship in place, the agreement between analytic 
theory and observations is remarkable. 

Heuristics. A fourth line of argument depends only upon the the mathematical structure 
of the TWE (in its undeveloped form relating S and Q) and the obeservation that the 
gradient of the SCZ angular velocity is dominated at mid-latitudes by its latitudinal angular 
component. We defer the details until the following section, when the TWE is explicitly 
presented. 

Any of the four arguments on its own would motivate a study of convection and differ- 
ential rotatation in which surfaces of constant angular velocity and constant residual entropy 
coincided. Taken as a whole, they comprise a truly compelling picture. 

We are thus motivated to extend the mathematical technique of B09 to fully convective 
stars, under the assumption that solar-like dynamical conditions prevail. In practice, this 
means solving the TWE not just for the 1/r gravitational potential appropriate to the SCZ 
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(r is the usual spherical radius), but for arbitrary potentials $(r). We present a suite of 
solutions for the resulting isorotation contours under a variety of simple assumptions for 
the mathematical form of the entropy/angular momentum relationship as well as for the 
type of boundary condition used for calculating Q. We consider interior angular velocity 
solutions for cases in which the surface rotation is monotonically increasing from pole to 
equator, monotonically decreasing in the same direction, and for zonal surface flows (seen, 
for example, in Jupiter and Saturn) as well. We also treat problems in which the angular 
velocity is given along the axis of rotation, in which case there may be very little surface 
variation. 

The profiles computed should provide a framework within which one may better un- 
derstand numerical simulations of rotation in convecting stars and planets. Moreover, if the 
rotation profiles of at least some fully convective stars are not dominated by magnetic fields, 
or turbulent or meridional fluxes, the calculated isorotation contours are viable predictions 
for the gross interior features of these bodies. Extracting such information from the global 
p-mode data of fully convective stars is a primary goal of the CoRoT satellite mission. 

An outline of our paper is as follows. Section 2 is a presentation of the mathematical 
solution for the isorotation contours of a fully convective star in which thermal wind balance 
is a good approximation. In section 3, we construct some specific solutions, corresponding to 
both solar-like and "antisolar" surface profiles, followed by models with banded zonal flows. 
Section 4 is a concluding discussion and summary. 



2. Analysis 

Throughout this paper, we follow the notation of BBLW. We denote cylindrical coor- 
dinates by radius R, azimuthal angle <fi and axial coordinate z. Spherical coordinates are 
given by (r, 9,<fr), where r is the radius from the origin, 9 is the colatitude angle, and is 
once again the azimuth. Unless otherwise stated, Q, the pressure P and the density p are 
understood to be azimuthal averages, independent of <fi. The dimensionless entropy function 
a is defined by: 

CT^lnPp- 7 , (1) 

where 7 is the usual (specific heat ratio) adiabatic index. As noted in the Introduction, 
$(r) will refer to the local gravitational potential. For the problem of interest, <3>(r) is a 
well-known, tabulated function. 

In thermal wind balance, one ignores contributions from convective turbulence (i.e., the 
convective Rossby number is small [Miesch & Toomre 2009]) and magnetic fields. This is an 
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important simplification that appears to work well for the bulk of the convective zone of the 
Sun (between roughly 0.75 and 0.95 solar radii), but must of course be reexamined critically 
for each application. 

The equation of thermal wind balance is (e.g. Thompson et al. 2003, Miesch 2005, 
B09): 

R ?f = (1) »£. (2) 

oz \jr ) dr o9 

For the SCZ, d<&/dr = GM & /r 2 , where G is the gravitational constant and M & represents a 
solar mass. As in BBLW, we assume that there is a functional relationship of the form 

a f = a-a r = f(n 2 ), (3) 

where a r is any function of r (in practice something very close to an angle-averaged a) and / 
is an unspecified function. We have noted earlier that this mathematical condition will have 
a physical basis if convection mixes the residual entropy a' in constant Q surfaces. But the 
mathematical structure of equation (j2D together with the most striking feature of the solar 
data already suggest a very simple heuristic argument (first mentioned in the Introduction) 
leading to a relation of the form seen in equation ([3]), an argument that is independent 
of whatever dynamical mechanism might be responsible for establishing the relationship. 
Helioseismology data show that at solar midlatitudes, the 9 gradient of Q dominates over its 
r gradient. Now in the TWE, a is operated upon only by the 9 gradient^. Thus, if a happens 
to have a dominant r gradient (not an unreasonable guess for the Sun), this component can 
largely be eliminated by subtracting off a suitably chosen radial function a r , i.e., by forming 
a a 1 quantity. Substituting a 1 for a does not change the TWE at all. But by its construction, 
this procedure will enhance the importance of the 9 gradient of a' relative to its r gradient, 
and therefore it motivates a search for solutions under the assumption that o' and Q share 
isosurfaces in common. This is precisely the content of equation ([3]) . 

It is of interest to consider why there can be such substantial departures from the 
Taylor-Proudman constraint of Q = Q(R) in a system that is so nearly adiabatic (and thus 
barotropic). The condition that isorotational and constant a' surfaces coincide provides an 
immediate answer: if O = Q(R) then a' = cr'(R) and da'/d9 must therefore be present. If 
this partial derivative is of order (RQ/v s ) 2 or larger, where v s is the adiabatic sound speed, 
it is enough to account for the order unity departure from Taylor-Proudman columns seen 
in the helioseismology data and our calculations. 



2 Likewise, Q appears only within a z gradient operator, and may have a function of R added to it without 
changing the equation. But this freedom is lost once £1 is specified on a particular surface. 



-6 - 



We assume that there is a stable, long lived average rotation profile Q(r,9) that is 
present in the star. The evolutionary processes that are responsible for the establishment of 
Q are likely to be complex, involving an as yet poorly understood combination of mechanical 
and thermal transport. For the problem at hand, however, these details are less important 
than the existence of Q(r, 9) itself, and the concomitant assumptions that the profile is in an 
average sense time steady and in thermal wind balance. If either of these assumptions fails 
for the fully convective stars we consider here, our model is no longer valid. 



Using equation (jBJ), the TWE O) becomes 



dQ 2 
Or 



f \ rf$ tanfl 



7r 2 sin 9 cos 9 J dr 



on 2 

i*r = ' (4) 



where /' = df /dQ 2 . The solution to equation @ is that Q 2 is constant on surfaces described 
by 

d9 ( f \ rf$ tan# 



This may be written 



dr \^fr 2 sin 9 cos 9 J dr 



(5) 



which has the solution 



— [r 2 sm 2 e + — $(r) = 0, (6) 
dr \ 1 J 

2 f 

r 2 sm 2 9 = A — $(r), (7) 

7 

where A is a constant of integration. This is a rather remarkable result, and worth stating 
more informally: if the TWE is obeyed and entropy is mixed as efficiently as possible in 
constant rotation surfaces, the isorotation contours are of the form R 2 = A + B<&(r), where 
A and B are constants along each contour. 

If our mathematical problem is to be well-posed, boundary conditions for Q 2 must 
be specified. One simple possibility is to demand that the angular velocity be a prescribed 
function of the colatitude angle 9q on some spherical surface r = ro- For the sake of simplicity, 
we consider this to be the surface of the star. Then the isorotation contour will be given by 

2f 

r 2 sin 2 9 = r 2 sin 2 9 - — [$(r) - $(r )]. (8) 

7 



In some stars, however, the characteristics may be quasi-spherical, in which case most of 
the isorotation contours will not penetrate the surface. (This seems to be true, for example, of 
convective stars with a poleward increasing surface rotation rate.) We require a formulation 
for the characteristics that is suitable for these cases. For quasi-spherical isorotation contours, 
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Q 2 is most conveniently prescribed as a function of zq along the axis of rotation. Then, the 
contours are defined by the equation 



r 2 sin 2 fl = — Mzo) - $(r)l. (9) 
7 

Since we will be concerned here with fully convective stars, the potential $(r) is that 
of an n = 1.5 polytrope, where n (the polytropic index) is given by 1/(7 — 1). To solve for 
$(r), recall that the equation of hydrostatic equilibrium may be immediately integrated for 
a polytropic star, giving (e.g. Schwarzschild 1958): 

7 P , , GM , . 

1 *M = , (10) 



7 - 1 p r 

where M is the mass of the star, and the integration constant on the right side is chosen 
so that the stellar potential matches the external potential —GM/r at the stellar surface 
r = vq. For the polytropic equation of state P = Kp 1 , where K is a constant defined by the 
adiabat, equation (TIP]) becomes 

(n + l)Kp 1/n + ®(r) = -—. (11) 

r 

To complete our determination of $(r), we need an explicit solution for p 1 /™ '. This is 
given by (Chandrasekhar 1939) 

p 1/n = pl /n Q n (0, (12) 
where p c is the central density of the star, £ is the dimensionless radius 

L a , = ^ + y'- (13 ) 

a 47rG 
and n (£) is the Lane-Emden function of order n, which satisfies 

d 2 (£6 



-e(e„) B . (14) 



Extensive tables of B n are available in the literature (e.g. Horedt 2004); for our needs a 
simple Pade approximant compiled by Pascual (1977) will suffice (see Appendix I). 

Returning to the problem at hand, solving equation (TXTT) for $ now gives 

$(r) = -— - (5/2)^e 3/2 (£). (15) 
^0 
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The zero of 63/2(0 is located at £0 = 3.65375 (Chandrasekhar 1939), corresponding to the 
surface r = r . Our final expression for the isorotation contour (jSJ) is given by 

r 2 sin 2 9 = rl sin 2 9 + 3fc 2 e 3/2 (£), (16) 

where c| (= KpJ n ) is the square of the isothermal sound speed at the center of the star, and 
we have used 7 = 5/3. The coefficient of O3/2 will be of order r\ if the dimensionless quantity 
Q 2 f (r f2/cs) 2 <C 1. The corresponding isorotation contour associated with equation 
is 

r 2 sin 2 fl = 3f 4[e 3/2 (0 - e 3/2 (zoWn))], (17) 

where z$ refers to the value of z at the start of the characteristic along the axis of rotation. 
The ratio z Q /r Q thus runs between and 1. 

For solar-like surface profiles, /' will generally be negative. This is because convection is 
more efficient along the rotation axis, so the poles will be associated with higher temperatures 
and small rotation rates relative to the equator. For antisolar profiles, similar reasoning 
suggests that /' should be positive. It is not yet possible to say more than this without 
detailed numerical simulations of the turbulent flow in convective stars. In the spirit of B09 
and BBLW, we therefore introduce the function F(sin 2 #o), which is of orer unity, rewriting 
equation (1TB]) as 

(r/r ) 2 sin 2 9 = (£/£ ) 2 sin 2 9 = sin 2 9 + F(sin 2 60)63/2(0, (18) 

where 

F(sm 2 9 ) = ^^<0 (19) 

is some simple dimensionless function of sin 2 6q. In practice, we will limit ourselves to linear or 
constant functions, there being at this stage no particular exigency for greater mathematical 
sophistication. 

For solar-like surface profiles, let us assume that F(sin 2 9 ) = — f3 2 sin 2 9 . Then, 
solving for sin 2 9 leads to 

sm e " = — 1 - WhM — ■ (20) 

If fl 2 is now a specified function at the surface, f2o(sin 2 9q) say, the solution throughout the 
stellar interior is given by substituting for sin 2 # from equation (I2U|) in the argument of Qq. 
This is the desired solution when Q is specified initially on a spherical surface. 



For the antisolar surface profiles, the zq inversion associated with equation ( I17p is some- 
what complicated, except when F is a (positive) constant, and we will restrict ourselves in 
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this work to this case. We then find 

(—) =A(X), X^e 3/2 (0-^sin^, (21) 

where the (very simple) A function, essentially an inversion of 63/2, is presented in Appendix 
I. Given the rotation profile along z , ^0(^0) sa y' the solution throughout the stellar interior 
is given by substituting for z\ from equation (}2"T]) in the argument of Q^. 

3. Contour description 

3.1. Solar surface profiles 

We choose the surface rotation profile f2(sin 2 9 ) to be equal to be that of the Sun, a 
simple Taylor expansion in sin 2 6q\ 

(tt /2n) = (386.2 + 198.7 sin 2 O - 66.7 sin 4 0„). (22) 

The units on the right are nano-Hz. Representative solutions are shown in figures 1-4. 

The four profiles of figures 1 and 2 correspond to a suite of increasing values of (3x, 
(0.2, 0.4, 0.6, 0.8) that show the contours evolving from nearly cylindrical Taylor-Proudman 
columns (corresponding to the limit fi\ — = 0), to a much more "flat-bottomed" con- 
figuration. In figures 3 and 4, we allow /3 X and (3 2 to vary separately, as indicated in the 
captions. Aside from a tightening or loosening of the contours, there is not a great qualitative 
difference between the cases of vanishing versus finite 02- 

The isorotation contours are, roughly speaking, divided into two classes, which we will 
term polar and equatorial. Polar contours extend from the surface, and cross through the 
rotation axis. Equatorial contours, by contrast, extend from the surface and cross through 
the equatorial plane. In some cases, there are contours that remain within the star, never 
reaching the surface. These tend to be quasi-spherical, deep in the interior near the stellar 
core. Mathematically, this corresponds to values of sin 2 9 in excess of unity, but the contours 
themselves are not unphysical: they simply remain below the surface of the star, and are 
valid solutions of the TWE. 

There is a particular value for 6$, say 8 C , which is the critical divide. Contours with 
#0 < 6 C are polar, while contours with 6 > 6 C are equatorial. In no cases do contours cross: 
caustics are absent from the solutions that we have studied. This is not true for the 1/r 
potential used in the calculations of B09 and BBLW. In that case, caustics formed when the 
solutions were extended beneath the convection zone. The Lane-Emden function $3/2 (0 5 
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Fig. 1. — Isorotation contours for a solar surface rotation profile. In all figures, the abscissa 
R and ordinate Z are given in units of the Lane-Emden radius a (eq. [13]). f3\ = 0.2, (3 2 = 
(left); ft = 0.4, 2 = 0.0 (right). 




Fig. 2.— As in figure(l), with f3 1 = 0.6, (3 2 = 0.0 (left); f3 1 = 0.8, /3 2 = 0.0 (right). 
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Fig. 4.— As in figure (1), with f3 1 = 0.5, f3 2 = 0.2 (left); ft = 0.2, j5 2 = 0.8 (right). 
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on the other hand, is evidently sufficiently well-behaved at small r to avoid the formation of 
caustics in the isorotation contours. 

Perhaps the most striking feature of these plots is how reminiscent of the Sun they are, 
despite the different gravitational potential and the extended domain. The contours show 
the familiar basic pattern of lying in planes of constant z near the axis of rotation, becoming 
much more cylindrical near the equator, and quasi-radial in the transition region between. 
Closer to the stellar surface, we see the characteristic solar feature of parallel lines once again 
in evidence. We have of course chosen a solar-like surface profile, but the results are generic 
for any slowly varying, poleward decreasing, surface profile with reflection symmetry about 
the equator. 



We next relax the constraint of a poleward decreasing surface rotation profile. Poleward 
increasing profiles are seen in simulations of convective red giant stars (Brun & Palacios 
2009). Since thermal convection is expected still to be most efficient parallel to the rotation 
axis (where it is unencumbered by the Coriolis force), the function f'(£l 2 ) is now negative. 
This causes interesting changes in the geometry of the isorotation contours. 



Figures (5) and (6) show a suite of internal isorotation contour diagrams for a poleward 
increasing surface profile: 



Compared with equation ( I2"2"j) for our fiducial solar poleward decreasing profile, we have 
changed the sign of the sin 2 9 term and reduced its magnitude by a factor of ten. We 
display four values of (3i, —0.5, —1, —1.5, —2, with f3 2 set equal to zero, so that /' is a global 
constant. The isorotation contours in this case do not resemble the Sun at all, for they 
now have a different topology. There is, for example, no critical contour dividing isorotation 
curves that start at the surface and then head either for the axis or for the equator. Instead, 
all contours cross the equatorial plane, and curve in the same sense (toward the rotation 
axis) as they rise upward from the equator. The most strongly inwardly curving contours 
first find the rotation axis (where they end), the less strongly curved contours encounter 
the surface before the axis. This is reminscent of the slowly rotating red giants in the Brun 
& Palacios (2009) simulations. In the numerical studies, however, the surface is in near 



3.2. Antisolar surface profiles. 



3.2.1. Vt initially specified on the surface 




(23) 
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Fig. 5. — Isorotation contours for an antisolar surface rotation profile (eq. [23]). As in figure 
(1), the abscissa R and ordinate Z are given in units of the Lane-Emden radius a (eq. p~3]). 
ft = -0.5, ft = (left); ft = -1.0, ft = 0.0 (right). 




Fig. 6. — As in figure (5), with ft = 
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uniform rotation and is best approached using the formalism of equation (Q (see below). 
For small values of the contours look, not surprisingly, cylindrical. This is clearly the 
Taylor-Proudman regime. But as the magnitude of fi\ increases, the contours become more 
spherical (prolate). In contrast to the case of solar surface profiles (and the Sun itself), in 
which the rotation rate increases with distance from the axis, the antisolar surface profiles 
are associated with an interior rotation structure with dfl 2 /dR < 0. 

3.2.2. Q initially specified along the rotation axis 

In stars showing significantly more Q variation along the axis of rotation than along a 
surface meridional contour, it makes sense to consider the initial rotation function to be of 
the form Q (zq), and use the formulation embodied in equations (fTTj) and (l2Tj) . We have 
adopted the model 

Q (z 2 ) = [l + ^oT 1 ' (24) 

leaving j3 as a parameter. Specifying the function F (treated here as a constant) and then 
uniquely defines the model. 

Figures (7) and (8) show a selection of results. In figure (7), we have set F = 2. On the 
left, /3 = 0.1; on the right, (5 = 1. In figure (8), F = 0.4 with (3 = 0.1 on the left and (3 = 1 
on the right. The larger value of F in figure (7) corresponds to a smaller rotation rate, and 
the isorotation contours are more rounded. Figure (8) corresponds to a larger rotation rate, 
and the contours are more elongated along the axis, distorted in the direction of Taylor- 
Proudman columns. Depending on the value of /3, the rotation can either be confined to 
a central core or extend to the surface. These findings compare very favorably with the 
morphologies found in the numerical simulations of Brun & Palacios (2009). 

3.3. Zonal surface flows 

Going further afield, one may speculate on what our theory would predict for the internal 
rotation profile of a convecting body whose surface exhibited zonal flow — as is seen in the 
giant planets of our solar system. This truly is a speculative domain, since even the sign 
of the da'/dQ 2 is not certain for this class of flows. Our motivation is simply to note an 
alternative to Taylor-Proudman columns for the internal flow contours of giant convecting 
gaseous planets and possibly brown dwarfs. In the case of Jupiter and Saturn the surface 
flows are likely to be shallow (e.g., Liu, Goldreich, & Stevenson 2008), and it is not known 
whether the convective interior supports such banded streams. Thus, whereas the solar 
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Fig. 7. — Internal isorotation contours for a convective star with specified along the 
rotation axis. The solution is generated by equations (ITTj) . ( 12T|) . and §M$- The abscissa R 
and ordinate Z are given in units of the Lane-Emden radius a (eq. [13]). F = 2 and /3 = 0.1 
(left) and 1 (right). 




Fig. 8.— As in figure (7), with F = 0.4, = 0.1 (left), (3 = 1 (right). 
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problem is tightly confined, and the deep convective envelope with a monotonic surface 
rotation profile is at least a well-motivated problem, this particular application is truly 
"untethered." These are purely formal solutions to the TWE. 

We consider a surface profile of the form 



In figures 9 and 10, we show a suite of four representative results with the equation (12D1 
parameters, fa = 0.2, 0.4, 0.6, 0.8 and fa = 0. As fa increases, the departure of the isorotation 
contours from Taylor- Proudman columns is evident, until by fa = 0.8 many contours show 
relatively little variation in y. Closer to the equator, however, the contours retain their 
cylindrical form. As noted in BBLW, rapid rotators are probably characterized by relatively 
small fa, because of the 1/fl 2 scaling of f'(Q 2 ). 

These results can be compared with numerical simulations of convectively driven dif- 
ferential rotation in Jupiter, Saturn, Uranus and Neptune. Such models have grown pro- 
gressively more realistic, moving from mean-field (Rieutord et al. 1984) to Boussinesq (e.g 
Aurnou & Heimpel 2004; Aurnou, Heimpel & Wicht 2007), and most recently to anelastic 
(Jones & Kuzanyan 2009) approximations. In the latter calculations, the Taylor-Proudman 
constraint is dominant, and the zonal velocity tends accordingly to be constant on cylinders. 



Our principal result, embodied in equation ([7]), is one of remarkable simplicity. If, 
in a convecting star, thermal wind balance holds and entropy is well mixed in surfaces of 
constant angular velocity, the form of the constant Q surface passing through a point in the 
r, 9 plane depends only upon the gravitational potential and a single "adjustable" constant. 
This constant is, in principle, determined from fundamental dynamics, but in practice it 
depends upon knowing how to treat the details of convective turbulent transport. One of 
the important consequences of this work, however, is precisely this point: ignorance of the 
details of turbulent transport does not mean complete ignorance of the consequences of the 
transport — in this case, the shapes of the isorotation surfaces. 

Unfortunately, unlike helioseismology, the fields of astero- and planetary seismology are 
still very much in their infancies. There are no data of quality comparable to the solar 
results with which we may compare our computed profiles. Until such data are available, 
the practical utility of our results will primarily be as an aid for understanding numerical 
simulations, and for theoretical investigations of convective processes in which a average 




(25) 



4. 



Discussion 
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Fig. 9. — Banded surface profile given by equation f l25|) . The abscissa R and ordinate Z 
are given in units of the Lane-Emden radius a (eq. [13]). 0\ = 0.2, 2 = (left); (3i = 0.4, 
02 = 0.0 (right). 




Fig. 10. — As in figure (9), with f3\ 



= 0.6, 02 = 0.0 (left); X = 0.8, 2 = 0.0 (right). 
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equilibrium rotation needs to be specified or understood (e.g. stellar evolution and dynamo 
theories). 

In principle, large scale numerical computations of fully convective stars are amenable 
to development at a level evinced by the solar convection studies (see, e.g. Miesch & Toomre 
2009). But in practice, this class of calculation is extremely demanding of resources, so that 
only the SCZ has been explored in detail and under a large variety of physical assumptions. 
Recent work, however, has shown that numerical simulations of convective envelopes are 
viable, and properly interpreted, a potentially rich source of dynamical information (e.g. 
Browning 2008; Brun & Palacios 2009 and references therein). The investigation of Brun 
& Palacios (2009) centered on four models of red giant stars. The more slowly rotating en- 
velopes produced quasi-spherical shells for the isorotational surfaces, while the more rapidly 
rotating models showed classic Taylor- Proudman cylinders. Our calculations also show this 
trend; indeed, strikingly so. More troublesome, for a direct application of the work pre- 
sented in the current paper, is the fact that our assumption of thermal wind balance was 
not a particularly accurate approximation to the vorticity equation in the simulations: both 
turbulent fluctuations and meridional circulation appeared to be more important in the red 
giant calculations than in the numerical simulations of the SCZ. Neither is convective Rossby 
number small for the simulations, as it is for the Sun. Yet the gross features of the red giant 
simulations seem rather well captured by our approach, and this is not because the models 
lack predictive power: inapproprately changing the sign of F, for example, destroys even the 
qualitative agreement with simulations in both the solar and red giant studies. Thus the 
marginality of our assumptions in the case of the latter may be more of a technical concern 
than a real conceptual impediment for understanding the qualitative structure of the flow. 
But clearly much more can be done. 

Another issue that we have not touched on is the effect of magnetic fields, shown in 
an extreme form by Browning's (2008) model of a fully convective M-star, where magnetic 
stresses transform a cylindrical f2-profile into almost uniform rotation. There is in fact an 
interesting analogy between the roles of residual entropy and of magnetic fields in coupling 
differentially rotating cylinders (Taylor 1963; Fearn 2007). 

Unlike the compelling fit of the solar isorotational contours presented in BBLW, the 
results presented in the current work must wait before observational comparisons are possible. 
Comparison with simulations is encouraging, however, and shows promise of genuine progress 
on both the analytic and numerical fronts. In it is particularly noteworthy, for example, 
that the poleward-increasing rotation profiles associated with at least some classes of fully 
convective stars show concave internal isorotation contours in simulations. We are in the 
early days yet of a field that is both rapidly and fruitfully developing. Even if only a small 
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class of fully convective systems has isorotation contours described by equation ([7j), this 
would be extremely useful. Not only would a piece of stellar hydrodynamics be secured, but 
a point of departure for more complex conditions would be established. 
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Appendix I: Pade approximants for © 3 / 2 (£)• 

For the calculations presented in this paper, we have adopted the following Pade ap- 
proximant for the n = 3/2 Lane-Emden function ©3/2(0 (Pascual 1977): 

03/2(0 * 1 + 6^ + ^ ' (26) 

where 

ax = 8.46561 x 10~ 2 , a 2 = 8.156966 x 10" 4 , (27) 

and 

61 = 8.201058 x 10" 2 , b 2 = 1.984127 x 10~ 3 . (28) 
A simpler, and only slightly less accurate approximation is given by 

©3/2(0 - (29) 

where 

a = 7.490699 x 10~ 2 , b = 9.175968 x 10~ 2 . (30) 



This preprint was prepared with the A AS IATgX macros v5.2. 
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The a and 6 constants are chosen so that the rational function goes to zero precisely at 
£ = £o = 3.65375, and to 1 — £ 2 /6 for small £. Contour plots produced by this approximation 
are almost indistinguishable from those shown in the paper. An advantage of this form for 
©3/2(0 is that it is easily invertible: 

e = = A(0), (31) 

a desirable quality in some applications. (We have suppressed the subscript 3/2 on 9 for 
clarity of presentation.) This defines the A function used in equation (1211) . 



Appendix II: Residual entropy bound 

The mathematical identity 



fda'\ ^f_dO_\ flda r 



\dr J g \d\nr J , \r 36 

may be used to establish a lower limit to the radial entropy gradient in the SCZ. Following 
the notion that constant a' and constant Q surfaces coincide, the first partial derivative on 
the right side of the equation may be read off directly from the helioseismology data. The 
second partial derivative on the right is evaluated via the TWE. The result is 

Da' \ n n fd\nQ\ ( 89 \ frQ 2 \ 

-2 7 tan — Kl— » ( 33 ) 



d In r J e \dlriz J R \dlar ) , \ g 

where g is the gravitational field of the SCZ, taken here to be 2.74 x 10 4 cm 2 s -1 . Using recent 
GONG data (kindly supplied to us by R. Howe), we estimate the logarithmic Q derivative 
as —0.12 and the 6-\nr gradient as —0.3. Then with r = 5.95 x 10 10 cm (or 0.85 solar radii), 
f2 = 2.5 x 10~ 6 s _1 , and 7 = 1.67, we obtain 

( 1 = -1.6 x 10" 6 tanfl. (34) 

\ a In r J e 

The numerical values correspond to regions of the SCZ near tan# = 1. Since a — a' + ay, 
this number is likely to be a lower bound (in magnitude) to the true midlatitude solar radial 
entropy gradient da/d In r. Mixing length theory gives values between 10 -6 and 10~ 5 for the 
absolute value of this gradient (e.g. Stix 2004). 



